{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flow and transport in 3D fractured medium\n",
    "\n",
    "Description: This notebook presents a weakly coupled flow and transport problem reported in the paper *PorePy: An Open-Source Simulation Tool for Flow and Transport in Deformable Fractured Rocks*, by\n",
    "Eirik Keilegavlen, Alessio Fumagalli, Runar Berge, Ivar Stefansson, Inga Berre. See arXiv:1712:00460 for details. The code in the notebook was that used to produce figure 5 in the paper, and if ran on a separate system, (almost) the same results should result - factors such as grid generation, linear solvers etc. may produce minor differences.\n",
    "\n",
    "The notebook doubles as an introduciton to PorePy for multi-physics, and is therefore rather long. To navigate quickly through the procedure, the main steps are: i) Create a mesh, ii) assign parameters for both flow and transport problems, iii) set up appropriate solvers, iv) discretize and solve.\n",
    "\n",
    "## Preliminaries\n",
    "To run this, you need to have PorePy installed and set up with Gmsh. Depending on the computer, the simulations can be somewhat memory and time consuming. To run on a subset of the fractures, use parameters for the grid generation function below.\n",
    "\n",
    "Also, it is recommended to install pyamg to speed up preconditioning for the pressure solver.\n",
    "\n",
    "## Imports\n",
    "The first step is to import all dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/Dropbox/workspace/python/porepy/src/porepy/grids/partition.py:18: UserWarning: Could not import pymetis. Some functions will not work as    intended\n",
      "  intended')\n"
     ]
    }
   ],
   "source": [
    "# External modules\n",
    "import numpy as np\n",
    "import pickle\n",
    "import scipy.sparse as sps\n",
    "import logging\n",
    "import sys\n",
    "\n",
    "# For plotting\n",
    "from IPython.display import HTML, display\n",
    "\n",
    "# import porepy\n",
    "import porepy as pp\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logger\n",
    "Next, set up a logger to monitor the progress.\n",
    "\n",
    "The syntax here is a bit confusing: The levels of the logging module in Python has names indicating that the program may break down (DEBUG, INFO, WARNING, ERROR, CRITICAL). We use them instead to set the level of verbosity, with DEBUG corresponding to 'print everything', and CRITICAL being almost silent. The tuning of these levels, and achieving consistency between modules, is still a work in progress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define\n",
    "root = logging.getLogger()\n",
    "root.setLevel(logging.WARN)\n",
    "\n",
    "# Only add handler if none is assigned. Without this, a new handler will be added every time this cell is ran.\n",
    "if not root.hasHandlers():\n",
    "    # Print to sys.stdout\n",
    "    ch = logging.StreamHandler(sys.stdout)\n",
    "    ch.setLevel(logging.INFO)\n",
    "    formatter = logging.Formatter('%(asctime)s - %(name)s - %(message)s')\n",
    "    ch.setFormatter(formatter)\n",
    "\n",
    "    root.addHandler(ch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grid generation\n",
    "The below function creates a 3D fracture etwork from a 2D outcrop. The extrusion produc vertical fractures only - more advanced options are available, but we have no geological basis for using them here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_grid(num_fractures=None):\n",
    "    \"\"\" Create a grid from an outcrop. \n",
    "    \n",
    "    To adjust the number of fractures to include, set num_fractures < 69 (maximum)\n",
    "    \"\"\"\n",
    "    # Default is to create all fractures\n",
    "    if num_fractures == None:\n",
    "        num_fractures = 66\n",
    "\n",
    "    # Fix the random seed, so that the produced network is reproducible\n",
    "    np.random.seed(7)\n",
    "    # Read points in outcrop from line, together with their connections.\n",
    "    # The file 'barely-connected_4pp.csv' is stochastically generated based \n",
    "    # on data from a real outcrop. Noteworthy, it contains several\n",
    "    # Y-intersections, which are preserved in extrusion.\n",
    "    # For more details, see tutorial on fracture extrusion.\n",
    "    pt, edges = pp.importer.lines_from_csv('barely-connected_4pp.csv')\n",
    "    \n",
    "    # Extrude fractures.\n",
    "    # Assign same family name to all fractures, and let them have zero incline.\n",
    "    # This is the place to modify if you want fractures with different \n",
    "    # properties, e.g. if they have different geological origins.\n",
    "    frac_family = np.zeros(edges.shape[1], dtype=np.int)\n",
    "    fractures = pp.extrusion.fractures_from_outcrop(pt, edges, family=frac_family,\n",
    "                                                    family_std_incline=[0.], tol=1e-5)\n",
    "\n",
    "    # Pick out the subset of fractures (or all) to be included in the grid\n",
    "    i = np.arange(num_fractures)\n",
    "    f = [fractures[fi] for fi in i]\n",
    "    # Create fracture network representation of the grid\n",
    "    network = pp.FractureNetwork(f)\n",
    "    # Assign tolerance used for computations. This is sometimes a critical \n",
    "    # parameter, it roughly defines how close to geometric objcets can be\n",
    "    # before they are considered identical (in case of two points), crossing\n",
    "    # (in case of two lines) etc.\n",
    "    # For general networks, it may be necessary to experiment a bit with this\n",
    "    # parameter to generate a reasonable mesh\n",
    "    network.tol = 1e-5\n",
    "    \n",
    "    # Export the network itself (no grid yet) as a vtu file, accesible by Paraview\n",
    "    network.to_vtk('network.vtu')\n",
    "\n",
    "    # Generate the mixed-dimensional grid.\n",
    "    # The mesh size is controlled by two parameters: mesh_size_frac is the target \n",
    "    # mesh size, while mesh_size_min is the minimal mesh size specified in the Gmsh\n",
    "    # configuration file. To which degree Gmsh adhers to this parameter\n",
    "    # varies from one case to the next.\n",
    "    # The parameter ensure_matching_face_cells=False is used to avoid an error\n",
    "    # message when matching cells with faces in the various dimensions\n",
    "    # The origin of the trouble seems to be in Gmsh (has been reported).\n",
    "    # We have never experienced any trouble with using this workaround, but\n",
    "    # some caution is still due.\n",
    "    gb = meshing.simplex_grid(network=network, mesh_size_frac=0.1, mesh_size_min=0.01,\n",
    "                              verbose=1,ensure_matching_face_cell=False)\n",
    "    \n",
    "    # gb is now a GridBucket, that is, a mixed-dimensional grid.\n",
    "    # Technically, this is implemented as a graph, with each node representing\n",
    "    # a grid (the matrix grid, a fracture, the intersection between fractures).\n",
    "    \n",
    "    # Rescale the domain from a meter to hundred meter scale, just to get more\n",
    "    # realistic numbers.\n",
    "    for g, _ in gb:\n",
    "        g.nodes *= 100 * pp.METER\n",
    "    gb.compute_geometry()\n",
    "    \n",
    "    return gb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem setup: Parameter specification and solvers\n",
    "To set up the flow and transport problems, we will use, respectively, an EllipticModel and a ParabolicModel. These are solvers dedicated to solving elliptic equations, e.g. the incompressible pressure equation, and advection-diffusion problems. The models automatically sets up simple discretization schemes, and tries to design decent linear solvers for the resulting systems of equations.\n",
    "\n",
    "### Pressure solver\n",
    "The first step is to provide simulation data. Simulation parameters are stored as part of the GridBucket, but accessing this can be somewhat cumbersome. To assist the assignment, and also to provide a reasonable way of setting default parameters, each of the models (pre-defined solvers), are accompanied by a DataAssigner. This can be used directly to define a simulaiton with default parameters, or modified as desired. Below, we give an example for the pressure equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The units module contains various physical constants\n",
    "# Set a relative high matrix permeability, this would correspond to\n",
    "# a quite high density of upscaled fractures.\n",
    "matrix_perm = 5 * pp.DARCY\n",
    "\n",
    "class MatrixDomainPressure(pp.EllipticDataAssigner):\n",
    "    \"\"\" Set data for the 3D domain (matrix) in the pressure equation.\n",
    "    \n",
    "    Fields that are not assigned here, will have the default values \n",
    "    prescribed in EllipticData (which again may point further to defaults\n",
    "    in the Parameter class).\n",
    "    \"\"\"\n",
    "    def __init__(self, g, dat):\n",
    "        self.dim = g.dim\n",
    "        pp.EllipticDataAssigner.__init__(self, g, dat)\n",
    "\n",
    "    # The permeability is set as isotropic and homogeneous.\n",
    "    # Both anisotropic and heterogeneous fields are easily accomodated\n",
    "    # by changing this line.\n",
    "    def permeability(self):\n",
    "        kxx = matrix_perm * np.ones(self.grid().num_cells)\n",
    "        return pp.SecondOrderTensor(self.dim, kxx)\n",
    "\n",
    "    # Boundary conditions: By default, Neumann conditions are assigned.\n",
    "    # Overwride these for two or the faces to set up flow along the y-axis.\n",
    "    # Note that this function only sets the type of boundary conditions,\n",
    "    # the exact values are set in the next function.\n",
    "    def bc(self):\n",
    "        g = self.grid()\n",
    "        # This will pick out faces with the maximum or minimum y coordinate\n",
    "        # We benefit from knowing that the simulation domain is a box alligned\n",
    "        # with the axes, if not, picking out faces would have been more cumbersome.\n",
    "        p_face = pp.params.bc.face_on_side(g, ['ymin', 'ymax'])\n",
    "        p_face = np.hstack((p_face[0], p_face[1]))\n",
    "        dir_str = p_face.size * ['dir']\n",
    "        return pp.BoundaryCondition(g, faces=p_face, cond=dir_str)\n",
    "\n",
    "    # Boundary values are set by specifying the value per face.\n",
    "    # If none are specified, homogeneous conditions are assumed.\n",
    "    # The reason for splitting types and values for boundary conditions is that \n",
    "    # types are needed at the time of discretization, while values are not imposed\n",
    "    # until assembley of linear system. \n",
    "    def bc_val(self):\n",
    "        g = self.grid()\n",
    "        y_min = pp.params.bc.face_on_side(g, 'ymin')\n",
    "        y_max = pp.params.bc.face_on_side(g, 'ymax')\n",
    "\n",
    "        bnd = np.zeros(g.num_faces)\n",
    "        bnd[y_min] = 0\n",
    "        bnd[y_max] = 50 * pp.BAR\n",
    "        return bnd\n",
    "\n",
    "    def porosity(self):\n",
    "        return 0.01 * np.ones(self.grid().num_cells)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, define parameters for the fractures, and their intersections. Most of the values are inherited from the matrix domain, but some are overwritten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FractureDomainPressure(MatrixDomainPressure):\n",
    "    def __init__(self, g, dat):\n",
    "        MatrixDomainPressure.__init__(self, g, dat)\n",
    "\n",
    "    # Define an aperture of 1cm.\n",
    "    # Note that the power is set according to the dimension, thus an \n",
    "    # intersection will have a transverse area of 1 cm squared.\n",
    "    def aperture(self):\n",
    "        return np.power(1 * pp.CENTI * pp.METER, 3 - self.grid().dim)\n",
    "\n",
    "    # Set a permeability of 10^-4, corresponding to the square of the aperture.\n",
    "    # The flow will again be scaled with the aperture inside the discretization,\n",
    "    # thus we in effect apply the so-called cubic law.\n",
    "    def permeability(self):\n",
    "        kxx = np.ones(self.grid().num_cells)*np.power(1e-4, self.grid().dim<3)\n",
    "        return pp.SecondOrderTensor(3, kxx)\n",
    "\n",
    "    # A lot can be said about the boundary conditions at fracture tips.\n",
    "    # We simply set a Neumann condition. Since this is the default, we \n",
    "    # need not define anything, but if we had not created the function, we would\n",
    "    # have inherited from the matrix domain. This might have been okay, as the\n",
    "    # fractures do not extend to the domain boundary, but the present version\n",
    "    # is cleaner and safer.\n",
    "    def bc(self):\n",
    "        return pp.BoundaryCondition(self.grid())\n",
    "\n",
    "    # No-flow conditions at fracture tips\n",
    "    def bc_val(self):\n",
    "        return np.zeros(self.grid().num_faces)\n",
    "\n",
    "    def porosity(self):\n",
    "        return np.ones(self.grid().num_cells)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having defined parameter classes for all geometric objects, assigning the data is easy: Simply loop over the GridBucket, and choose DataAssigner according to the grid dimension.\n",
    "\n",
    "If the data was heterogeneous, for instance some fractures are blocking, this would have required a few more steps. First, an extra DataAssigner is needed (BlockingFractureDomainPressure?), which would inherit from FractureDomainPressure, and likely override the permeability() function. Second, we need a way to differ between the individual fracture objects. Fracture grids will have an attribute frac_num, which reflects the order of the fractures when they were fed into the FractureNetwork, and further into the meshing algorithm. This can be used to identify the individual fractures.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define method to assign parameters to all nodes in the GridBucket\n",
    "def assign_pressure_data(gb):\n",
    "    # Prepare the model to recieve problems to be solved.\n",
    "    gb.add_node_props(['problem'])\n",
    "    # Loop over the grid bucket. This loops over the grid, and the\n",
    "    # data field associated with the grid. The latter can be used to stored\n",
    "    # pretty much anything, among them the problem definition.\n",
    "    for g, d in gb:\n",
    "        # The 3D domain is assigned Matrix properties\n",
    "        if g.dim == 3:\n",
    "            d['problem'] = MatrixDomainPressure(g, d)\n",
    "        # All others get fracture properties.\n",
    "        elif g.dim < 3:\n",
    "            d['problem'] = FractureDomainPressure(g, d)\n",
    "    # Assign coupling discharge\n",
    "    gb.add_edge_props(['param'])\n",
    "    for e, d in gb.edges():\n",
    "        g_h = gb.nodes_of_edge(e)[1]\n",
    "        d['param'] = pp.Parameters(g_h)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, define the pressure solver. It would have been equally natural to let assign_pressure_data() be incorporated into __init__, but we pulled it out here to give it the proper attention:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PressureSolver(pp.EllipticModel):\n",
    "    \"\"\" Solver for the pressure equation.\n",
    "    \n",
    "    No methods are specified here, thus all functionality is inherited from \n",
    "    EllipticModel. \n",
    "    \"\"\"\n",
    "    def __init__(self, gb, solver_data):\n",
    "        assign_pressure_data(gb)\n",
    "        pp.EllipticModel.__init__(self, gb)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transport discretization\n",
    "So far, we have defined a solver for the pressure equation only. Next, do the same for the transport solver. The steps are the same: To assign parameters, we will rely upon the ParabolicDataAssigner for default values, and overwrite values when appropriate. Then define a solver, and we are done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MatrixDomainTransport(pp.ParabolicDataAssigner):\n",
    "    def __init__(self, g, dat):\n",
    "        self.dim = g.dim\n",
    "        pp.ParabolicDataAssigner.__init__(self, g, dat)\n",
    "\n",
    "    # Thermal diffusivity of granite\n",
    "    def diffusivity(self):\n",
    "        kxx = 3 * np.ones(self.grid().num_cells)\n",
    "        return pp.SecondOrderTensor(self.grid().dim, kxx)\n",
    "\n",
    "    # Boundary conditions are set similar to the pressure equation:\n",
    "    # Impose a gradient in the y-direction, and insulate the other boundaries.\n",
    "    def bc(self):\n",
    "        g = self.grid()\n",
    "        p_face = pp.params.bc.face_on_side(g, ['ymin', 'ymax'])\n",
    "        p_face = np.hstack((p_face[0], p_face[1]))\n",
    "        dir_str = p_face.size * ['dir']\n",
    "        return pp.BoundaryCondition(g, faces=p_face, cond=dir_str)\n",
    "\n",
    "    # Boundary values are set to 0 and 1. Can easily be changed by scaling.\n",
    "    def bc_val(self, t):\n",
    "        g = self.grid()\n",
    "        y_min = pp.params.bc.face_on_side(g, 'ymin')\n",
    "        y_max = pp.params.bc.face_on_side(g, 'ymax')\n",
    "\n",
    "        bnd = np.zeros(g.num_faces)\n",
    "        bnd[y_min] = 0\n",
    "        bnd[y_max] = 1\n",
    "        return bnd\n",
    "\n",
    "    # For implementation reasons, domains of all dimensions have an apperture.\n",
    "    # In fractures and fracture intersections, this is used to scale cell volumes to a \n",
    "    # Nd value. In 3D, simply set the aperture to 1; this corresponds to a unit scaling.\n",
    "    def aperture(self):\n",
    "        return np.ones(self.grid().num_cells)\n",
    "\n",
    "    def porosity(self):\n",
    "        return 0.01 * np.ones(self.grid().num_cells)\n",
    "\n",
    "    # Rock density corresponding (approximate) average values for granite. Values from Wikipedia.\n",
    "    def rock_density(self):\n",
    "        return 2.7 * pp.KILO * pp.KILOGRAM / pp.METER**3\n",
    "\n",
    "    # Rough estiamet of water density\n",
    "    def fluid_density(self):\n",
    "        return 1.0 * pp.KILO * pp.KILOGRAM / pp.METER**3\n",
    "\n",
    "    # Rock heat capasity for granite\n",
    "    def rock_specific_heat(self):\n",
    "        return 790. / pp.KILOGRAM\n",
    "\n",
    "    # Heat capasity corresponds to water in standard conditions\n",
    "    def fluid_specific_heat(self):\n",
    "        return 4185. / (pp.KILOGRAM)\n",
    "\n",
    "\n",
    "class FractureDomainTransport(MatrixDomainTransport):\n",
    "    \"\"\" The FractureDomain inherits from Matrix domain. The difference is in the \n",
    "    aperture, that has the same value as for the flow problem.\n",
    "    \n",
    "    \n",
    "    \"\"\"\n",
    "    def __init__(self, g, dat):\n",
    "        MatrixDomainTransport.__init__(self, g, dat)\n",
    "\n",
    "    def aperture(self):\n",
    "        return np.power(1 * pp.CENTI * pp.METER, 3 - self.grid().dim)\n",
    "\n",
    "    def bc(self):\n",
    "        # Overwrite MatrixDomain's bc, we're using only Neumann\n",
    "        return pp.params.bc.BoundaryCondition(self.grid())\n",
    "\n",
    "    def bc_val(self, t):\n",
    "        return np.zeros(self.grid().num_faces)\n",
    "\n",
    "    def porosity(self):\n",
    "        return np.ones(self.grid().num_cells)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, define the transport solver. The most noteworthy point here is the solver() method, which returns an implicit time stepper with a linear solver tailored for large, advection-dominated problems. Also, this is the place where the end time for simulation and time step size is set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TransportSolver(pp.ParabolicModel):\n",
    "    def __init__(self, gb):\n",
    "        \"\"\" Here the data assignment is located in __init__,\n",
    "        as indicated in the pressure solver\n",
    "        \"\"\"\n",
    "        gb.add_node_props(['transport_data'])\n",
    "        for g, d in gb:\n",
    "            if g.dim == 3:\n",
    "                d['transport_data'] = MatrixDomainTransport(g, d)\n",
    "            elif g.dim < 3:\n",
    "                d['transport_data'] = FractureDomainTransport(g, d)\n",
    "\n",
    "        pp.ParabolicModel.__init__(self, gb)\n",
    "\n",
    "    # The spatial discretization returns both advective and diffusive\n",
    "    # discretization. If we were to solve a pure advective problem, we\n",
    "    # would simply return only self.advective_disc()\n",
    "    def space_disc(self):\n",
    "        return self.advective_disc(), self.diffusive_disc()\n",
    "    \n",
    "    def solver(self):\n",
    "        # Create an ad hoc implicit time stepper that inherits from the\n",
    "        # standard Implicit class (backward Euler approach).\n",
    "        # The only change we do is to let the linear solver use\n",
    "        # GMRES, preconditioned with ILU for large problems.\n",
    "        # The definition of 'Large' is somewhat arbitrary; the point is \n",
    "        # to turn to an iteratie solver for systems large enough that a\n",
    "        # direct solver is inefficient, or consumes too much memory.\n",
    "        # What that means is very much a user specific parameter; \n",
    "        # the value prescribed here should be conservative in most cases.\n",
    "        class ImplicitSolverGmres(pp.Implicit):\n",
    "\n",
    "            def __init__(self, model):\n",
    "                pp.Implicit.__init__(self, model)\n",
    "\n",
    "            def step(self):\n",
    "                from porepy.numerics.linalg.linsolve import Factory as LSFactory\n",
    "                A = self.lhs\n",
    "                b = self.rhs\n",
    "\n",
    "                ls = LSFactory()\n",
    "                if b.size < 40000:\n",
    "                    self.p = ls.direct(A, b)\n",
    "                    return self.p\n",
    "                else:\n",
    "                    # Note use of ILU here: This is appropriate for \n",
    "                    # advection dominated problems, where most of the action\n",
    "                    # is local. For conduction dominated problems, a better\n",
    "                    # approach would be to use elliptic solvers, e.g. a\n",
    "                    # algebraic multigrid approach\n",
    "                    precond = ls.ilu(A)\n",
    "                    slv = ls.gmres(self.lhs)\n",
    "                    # The use of callback here will log information on the\n",
    "                    # iteration procedure.\n",
    "                    # Parameters for maxiter, restart etc. are somewhat random.\n",
    "                    # In practice, the iterations converges quite fast for this \n",
    "                    # problem.\n",
    "                    self.p, info = slv(self.rhs, M=precond, callback=True,\n",
    "                               maxiter=10000, restart=1500, tol=1e-8)\n",
    "                    return self.p\n",
    "        return ImplicitSolverGmres(self)\n",
    "\n",
    "    def end_time(self):\n",
    "        return 10 * pp.YEAR\n",
    "\n",
    "    def time_step(self):\n",
    "        return self.end_time()/10\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up solvers, then solve\n",
    "We are finally ready to define our solver objects and solve for flow and temperature. With all parameters defined, this is a relatively simple code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/anaconda3/envs/porepy/lib/python3.5/site-packages/vtk/util/numpy_support.py:134: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n",
      "  assert not numpy.issubdtype(z.dtype, complex), \\\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Construct mesh\n",
      "Use existing decomposition\n",
      "Minimal distance between points encountered is 0.004273543932074861\n",
      "Gmsh processed file successfully\n",
      "\n",
      "\n",
      "Grid creation completed. Elapsed time 0.1820986270904541\n",
      "\n",
      "\n",
      "Created 1 3-d grids with 10989 cells\n",
      "Created 3 2-d grids with 1100 cells\n",
      "\n",
      "\n",
      "Done. Elapsed time 1.0201530456542969\n",
      "Assemble in bucket\n",
      "Done. Elapsed time 0.026569604873657227\n",
      "Compute geometry\n",
      "Done. Elapsed time 0.047032833099365234\n",
      "Split fractures\n",
      "Done. Elapsed time 0.7815241813659668\n",
      "Mesh construction completed. Total time 1.8778269290924072\n"
     ]
    }
   ],
   "source": [
    "# Define computaitonal mesh (mixed-dimensional)\n",
    "# You may also get a warning here, related to inconsistency between cells and faces.\n",
    "# This is caused by what appears to be a bug in gmsh (we have reported),\n",
    "# that shows up for small and skewed cells and faces.\n",
    "# We have not experienced any practical implications of this.\n",
    "mesh = create_grid()\n",
    "\n",
    "# Create an exporter object, and dump the grid\n",
    "# NOTE: The export, and particular representing the grid in a format \n",
    "# amenable for the unstructured data format in vtu is time consuming\n",
    "# for a pure python code. If numba is available (conda install numba),\n",
    "# a just in time compiler will be invoked, and the whole process\n",
    "# speeds up dramatically. Do install.\n",
    "mesh_exp = pp.Exporter(mesh, 'grid')\n",
    "mesh_exp.write_vtk()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting grid looks like this, after some manipulation in Paraview"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><td><img src='fig/g0.png'></td><td><img src='fig/g1.png'></td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML(\"<table><tr><td><img src='fig/g0.png'></td><td><img src='fig/g1.png'></td></tr></table>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pressure solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2018-03-19 14:37:21,083 - porepy.numerics.elliptic - Solve elliptic model\n",
      "2018-03-19 14:37:21,084 - porepy.numerics.elliptic - Discretize\n",
      "2018-03-19 14:37:21,165 - porepy.numerics.elliptic - Done. Elapsed time 0.08112049102783203\n",
      "2018-03-19 14:37:21,165 - porepy.numerics.elliptic - Solve linear system using direct solver\n",
      "2018-03-19 14:37:21,352 - porepy.numerics.elliptic - Done. Elapsed time 0.18686413764953613\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/anaconda3/envs/porepy/lib/python3.5/site-packages/vtk/util/numpy_support.py:134: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n",
      "  assert not numpy.issubdtype(z.dtype, complex), \\\n"
     ]
    }
   ],
   "source": [
    "# Define pressure solver for the given grid.\n",
    "# This will assign parameters, using the above classes\n",
    "pressure_solver = PressureSolver(mesh, None)\n",
    "\n",
    "# Solve the pressure equation.\n",
    "#\n",
    "# TECHNICAL INFORMATION:\n",
    "# This discretizes the equation using a two-point flux approximation \n",
    "# (unless the discretization method is overwritten), assembles \n",
    "# the linear system and right hand side contributions from sources\n",
    "# and boundary conditions. The linear system is then solved, hopefully\n",
    "# with a somewhat efficient solver. Exactly how this should be defined is \n",
    "# not entirely clear at the moment (the linear system has saddle-point\n",
    "# characteristics, which makes solution a bit challenging). What we do\n",
    "# for now is to apply a direct solver for small systems (specified by the\n",
    "# parameter max_direct). For large problem, a domain decomposition approach\n",
    "# is applied for the rock matrix and all fractures and intersection. For\n",
    "# large sub-problems (in practice the matrix), an AMG solver is invoked\n",
    "# if available, if not, ILU is the fallback.\n",
    "# It is strongly recommended to install pyamg before solving large problems.\n",
    "#\n",
    "# UPPSHOT:\n",
    "# For large problems, this may take some time, and you may want to monitor memory usage.\n",
    "\n",
    "pressure_solver.solve(callback=True)  # Use callback to monitor progression\n",
    "\n",
    "# The pressure solution is stored as a single vector, considering all\n",
    "# geometric objects (rock matrix, fractures, intersections) at the same time\n",
    "# Split into the different components in the computational mesh.\n",
    "pressure_solver.split()\n",
    "\n",
    "# Store the solution in the data-part of the computational grid.\n",
    "# This saves the solution to vtk.\n",
    "pressure_solver.save('solution')\n",
    "\n",
    "# Derive fluxes from the computed pressures\n",
    "pp.fvutils.compute_discharges(mesh, p_name='solution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is what the pressure solution looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table><tr><td><img src='fig/pressure.png'></td><td><img src='fig/p2.png'></td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML(\"<table><tr><td><img src='fig/pressure.png'></td><td><img src='fig/p2.png'></td></tr></table>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transport\n",
    "And finaly, the transport solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2018-03-19 14:38:09,893 - porepy.numerics.time_stepper - Time stepping using 10 steps\n",
      "2018-03-19 14:38:09,977 - porepy.numerics.time_stepper - Step 0 out of 10\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eke001/anaconda3/envs/porepy/lib/python3.5/site-packages/vtk/util/numpy_support.py:134: FutureWarning: Conversion of the second argument of issubdtype from `complex` to `np.complexfloating` is deprecated. In future, it will be treated as `np.complex128 == np.dtype(complex).type`.\n",
      "  assert not numpy.issubdtype(z.dtype, complex), \\\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2018-03-19 14:38:10,406 - porepy.numerics.time_stepper - Step 1 out of 10\n",
      "2018-03-19 14:38:10,824 - porepy.numerics.time_stepper - Step 2 out of 10\n",
      "2018-03-19 14:38:11,371 - porepy.numerics.time_stepper - Step 3 out of 10\n",
      "2018-03-19 14:38:11,786 - porepy.numerics.time_stepper - Step 4 out of 10\n",
      "2018-03-19 14:38:12,202 - porepy.numerics.time_stepper - Step 5 out of 10\n",
      "2018-03-19 14:38:12,721 - porepy.numerics.time_stepper - Step 6 out of 10\n",
      "2018-03-19 14:38:13,147 - porepy.numerics.time_stepper - Step 7 out of 10\n",
      "2018-03-19 14:38:13,613 - porepy.numerics.time_stepper - Step 8 out of 10\n",
      "2018-03-19 14:38:14,106 - porepy.numerics.time_stepper - Step 9 out of 10\n"
     ]
    }
   ],
   "source": [
    "# Define transport solver\n",
    "transport_solver = TransportSolver(mesh)\n",
    "\n",
    "# Discretize and perform time stepping.\n",
    "# The time stepper discretizes the advective term using a \n",
    "# first order upwind approach, while the conductive term is treated\n",
    "# by TPFA.\n",
    "# We saw above how the linear solver was defined in this case.\n",
    "# To save each time step we send in the key-argument save_as, which\n",
    "# will save the solution in a file \"results/solution.pvd\" as the variable\n",
    "# \"transport\".\n",
    "sol = transport_solver.solve(save_as='transport')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we are done. Visualized, the transport solution evolves like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HTML('<img src=\"fig/out.gif\">')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comments and outlook\n",
    "The notebook turned out to be rather lengthy, but if we strip away all more or less helpful comments, the following steps remains:\n",
    "1. Define computational mesh.\n",
    "2. Define all parameters necessary for the solution. By default, these will typically be of unit size; if other values are needed they must be overwritten. PorePy also have a convenient way of differing between parameter assignment for the rock matrix and (individual) fractures.\n",
    "3. Implement a discretization strategy. PorePy defines standard solvers for pressure and transport problems, using simple but robust discretization techniques. The standard models can be modified when needed, typically in the case where the user knows something about the structure of the problem at hand.\n",
    "4. Solve the problems, save relevant data.\n",
    "\n",
    "Needless to say, these steps are part of any solution strategy.\n",
    "\n",
    "## Possible extensions\n",
    "As the intention of the example was to illustrate multi-physics couplings in the setting of flow and transport, we comment on two natural extensions:\n",
    "1. Account for fluid compressibility in the flow field. Assuming the fluid is water, this can be modeled by replacing the EllipticModel with a SlightlyCompressibleModel. As the flow field is now time dependent, it should be updated during the transport solve. Exactly what is the most efficient approach there depends on the coupling strength - operator splitting is a natural option. Anyhow, it would be necessary to explictly write out the time loop containing pressure and transport solves. Apart from this most of the script should be the same.\n",
    "2. As an extension of the first point, in certain applications (e.g. geothermal energy) in may be of interest to account for more ellaborate couplings, such as temperature dependency of the viscosity etc. This will in general introduce non-linear couplings that, if strong, should be solved simultaneously. This is beyond the current scope of PorePy, but should be feasible by exploiting free python libraries for automatic differentiation.\n",
    "\n",
    "Another, more technical, step is to assign heterogeneous parameter values, mostly for cell-values fields such as permeabilities. This would require functions to define parameters e.g. based on cell center coordinates. Once this is in place, it is simply a matter of overwriting the relevant functions in data assigners such as MatrixPressureDomain etc."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
